* What does this do? Create Figure 4 and Figure A9

clear all

use "$savedata/masterdata.dta", replace

* 3-year experience of non-STEMI cases
gen nstemi3 = ami3-stemi3

keep if sample25==1
gen vol = vol25

*Split into stemi and non-stemi cases
gen ps2 = 2
replace ps2 = 1 if stemi2==1 
table ps2, c(mean all_death30 mean all_death365)

foreach x in 1 2{
gen one`x' = 0
replace one`x' = 1 if ps2==`x'

egen n_efe_n`x' = sum(one`x'), by(doctor_id)
egen vol`x' = sum(one`x'), by(doctor_id)
}

foreach x in 1 2{
gen one`x'_2017 = 0 
replace one`x'_2017 = 1 if ps2==`x' & finyear==2017
}


foreach x in 1 2{
reghdfe survive365 c.prevyear_cost if ps2==`x', absorb(i.sex##i.derv_age i.black i.mixed i.chinese i.asian i.race_miss i.ynch* i.prevyear_stroke i.di1 i.di2 i.di3 i.di4 i.di5 i.shock i.arythmia i.arthero i.arrest i.dow##i.admidate_mont##i.finyear i.hyid, savefe) keepsingleton
predict residual`x', residuals

xtreg residual`x' c.stemi3 c.nstemi3 c.std_nonami3, fe i(doctor_id)
predict docfe_n`x', u
gen sigma_u_n`x' = e(sigma_u)
gen sigma_e_n`x' = e(sigma_e)
}

foreach x in 1 2{
gen signal_n`x' = ((sigma_u_n`x'^2) / ((sigma_u_n`x'^2) + ((sigma_e_n`x'^2)/n_efe_n`x')))
gen adj_docfe_n`x' = docfe_n`x'*signal_n`x'
}


* Collapse to doctor level
collapse (mean) adj*, by(doctor_id)

su adj* /* gives the standard deviation of quality for each type reported in Section 6.1 */

* Check correlation between two measures (reported in Section 6.1)
corr adj_docfe_n1 adj_docfe_n2

* Make Figure 4
twoway kdensity adj_docfe_n1, lcolor(black) lpattern(dash) || kdensity adj_docfe_n2, lpattern(solid) lcolor(black) ///
legend(label(1 "STEMI (S)") label (2 "Non-STEMI (N)") row(1)) ///
graphregion(color(white)) ytitle("Density") xtitle("Estimated patient type-specific doctor FE") xlab(-0.15(0.05)0.15)
graph export "$results/Figure4.pdf", as(pdf) replace

* Make Figure A9
* Compare scores for each type
scatter adj_docfe_n1 adj_docfe_n2, mcolor(gs3) || lfit adj_docfe_n1 adj_docfe_n2, legend(off) ytitle("Estimates of {&mu}{superscript:S}") xtitle("Estimates of {&mu}{superscript:N}") graphregion(color(white)) 
graph export "$results/FigureA9.pdf", as(pdf) replace

